home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / piw131.zip / PIW.TXT < prev    next >
Text File  |  1994-07-29  |  9KB  |  197 lines
  1.  
  2. PiW - Compute Pi to a million or so digits                  29 JULY 1994
  3.  
  4. The program PiW can compute Pi to a million or so decimal places. It is
  5. written in Borland C++ for Windows, so it requires an IBM PC compatible
  6. computer running Microsoft Windows 3.0 or later. It is recommended that
  7. Windows be run in the 386 enhanced mode. With 16 megabytes of RAM it can
  8. compute Pi to 524,200 decimal places, with 32 megabytes of RAM 1,048,500
  9. places. Run the executable file PiW.Exe and you will get the message:
  10.  
  11.   PiW - Compute Pi to a million or so decimal places in Windows
  12.   Version 1.31, last revised: 1994-07-29, 0600 hours
  13.   Copyright (c) 1981-1994 by author: Harry J. Smith,
  14.   19628 Via Monte Dr., Saratoga, CA 95070.  All rights reserved.
  15.  
  16.   Will write file: PiW.Out containing Pi.
  17.  
  18.   Will write file: PiW.Log containing a log of operations.
  19.  
  20.   Will write file: PiReset.Txt containing a restart file.
  21.  
  22.   Will write file: PiReset.Old containing a restart file backup.
  23.  
  24.   How many decimal places do you need? (n < 0 to restart) (0 to quit): _
  25.  
  26. Type in a number like 2000 and press Enter. The program will respond
  27. with:
  28.  
  29.   Use algorithm b? (Y/N): _
  30.  
  31. Type a single letter (y, Y, n, or N) and press Enter. The program will
  32. respond with:
  33.  
  34.   How many digits between commas in output? (0 for no , < 0 for digits only): _
  35.  
  36. Three would be the most common input here, but zero and five are also
  37. common. A negative input causes nothing but the digits of Pi to be
  38. output. An input of -32768 is treated as minus zero.
  39.  
  40. Type a number and press Enter. The program will start computing Pi and
  41. output a log of its operations with diagnostic and timing data to the
  42. screen and to the file PiW.Log. If you answer with a character other
  43. than y or Y above the program will use algorithm a instead of algorithm
  44. b. For a request of 2000 decimal places, algorithm b, and no commas, the
  45. output looks like:
  46.  
  47.   Start of PiW.Cpp 1994-07-29 06:30:21.49
  48.   Trying to compute Pi to about 2023 decimal places.
  49.   T = 0  DT = 0 sec.  Starting Pi
  50.     Starting to compute Pi by algorithm b
  51.   T = 2.31  DT = 2.31 sec.  Starting Iteration
  52.     Start of iteration number 1 of 5
  53.   T = 19.17  DT = 16.86 sec.  Starting Iteration
  54.     Start of iteration number 2 of 5
  55.   T = 34.94  DT = 15.77 sec.  Starting Iteration
  56.     Start of iteration number 3 of 5
  57.   T = 52.68  DT = 17.74 sec.  Starting Iteration
  58.     Start of iteration number 4 of 5
  59.   T = 70.64  DT = 17.96 sec.  Starting Iteration
  60.     Start of iteration number 5 of 5
  61.   T = 89.26  DT = 18.62 sec.  All done but...
  62.     All done but inverting alpha
  63.   T = 91.67  DT = 2.41 sec.  Pi is done
  64.     All done computing Pi
  65.  
  66.   T = 91.73  DT = 0.06 sec.  Writing Pi to disk
  67.  
  68.   Pi = m.n E+0, m.n =
  69.  
  70. 3.141592653589793238462643383279502884197169399375105820974944592307816406286208
  71. 99862803482534211706798214808651328230664709384460955058223172535940812848111745
  72. 02841027019385211055596446229489549303819644288109756659334461284756482337867831
  73. 65271201909145648566923460348610454326648213393607260249141273724587006606315588
  74. 17488152092096282925409171536436789259036001133053054882046652138414695194151160
  75. 94330572703657595919530921861173819326117931051185480744623799627495673518857527
  76. 24891227938183011949129833673362440656643086021394946395224737190702179860943702
  77. 77053921717629317675238467481846766940513200056812714526356082778577134275778960
  78. 9173637178721468440901224953430146549585371050792279689258923 E+0 (700) [2051]
  79.  
  80.   T = 93.82  DT = 2.09 sec.  All done
  81.   End of PiW.Cpp 1994-07-29 06:31:55.41
  82.  
  83. It is recommended that, if you compute Pi to more places than you can
  84. verify to be correct by an independent source, you should compute Pi to
  85. this number of places using both algorithm a and algorithm b and use the
  86. DOS program FC (File Compare) to verify that the same value was computed
  87. by both algorithms.
  88.  
  89. These algorithms are documented in Scientific American, February 1988,
  90. Ramanujan and Pi, by Jonathan M. Borwein and Peter B. Borwein.
  91.  
  92. Algorithm a:
  93.  
  94.   Let y[0] = SqRt(1/2),    x[0] = 1/2
  95.  
  96.   y[n] = (1 - SqRt(1 - y[n-1]^2)) / (1 + SqRt(1 - y[n-1]^2))
  97.  
  98.   x[n] = ((1 + y[n])^2 * x[n-1]) - 2^n * y[n]
  99.  
  100. Algorithm b:
  101.  
  102.   Let y[0] = SqRt(2) - 1,    x[0] = 6 - 4 * SqRt(2)
  103.  
  104.   y[n] = (1 - SqRt(SqRt(1 - y[n-1]^4))) / (1 + SqRt(SqRt(1 - y[n-1]^4)))
  105.  
  106.   x[n] = ((1 + y[n])^4 * x[n]) - 2^(2n+1) * y[n] * (1 + y[n] + y[n]^2)
  107.  
  108. For both algorithms, xn converges to 1/Pi. Algorithm a is quadratically
  109. convergent and algorithm b is quartically convergent. The following
  110. table shows how many iterations are needed to compute Pi to a given
  111. number of significant digits:
  112.  
  113.   Iterations  Iterations    Digits
  114.   of algo. a  of algo. b  matching Pi
  115.   ----------  ----------  -----------
  116.       1                        0
  117.       2                        3        (1/x2 = 3.140, algo. a)
  118.       3           1            8
  119.       4                       19
  120.       5           2           41
  121.       6                       84
  122.       7           3          171
  123.       8                      345
  124.       9           4          694
  125.      10                     1392
  126.      11           5         2788
  127.  
  128. Algorithm b converges to a given number of digits about twice as fast as
  129. algorithm a, but takes about twice as much work per iteration. Using
  130. PiW, algorithm b is 10 to 20 percent faster than algorithm a.
  131.  
  132. I used PiW to compute Pi to 500,000 decimal digits. This was done on an
  133. IBM AT compatible 33 MHz 486 DX computer using Windows 3.1 with 16
  134. megabytes of RAM and 22 mega bytes of virtual memory. It took 37.3 hours
  135. for algorithm a and 31.0 hours for algorithm b. The results were the
  136. same. The divide and square root routines have been improved since these
  137. runs, so the program is a little faster now (10 to 20%).
  138.  
  139. The distribution disk for PiW includes the executable file PiW.Exe, all
  140. the source code, this document, and the file PiW500K.Out containing Pi
  141. computed to 500,000 decimal digits. Read the READ-ME and the WHATFOR
  142. files for a complete description of what's on the disk.
  143.  
  144. The first time I tried to compute Pi to more places than is normally
  145. needed I computed it to 50 decimal places using paper, pencil, and a
  146. hand calculator. The equation used was
  147.  
  148.         Pi = 16 * ArcTan(1/5) - 4 * ArcTan(1/239)
  149.  
  150.         ArcTan(1/x) = 1/x - 1/(3 * x^3) + 1/(5 * x^5) - ...
  151.  
  152. This was done in the early 1970's. The equation was derived in 1706 by
  153. John Machin (1685-1751).  After I got my Apple II+ computer in 1979, I
  154. computed Pi to 15,300 decimal places using Apple Pascal and the equation
  155.  
  156.   Pi = 48 * ArcTan(1/18) + 32 * ArcTan(1/57) - 20 * ArcTan(1/239)
  157.  
  158. This is the same equation I used in 1989 to compute a 114,632 decimal
  159. place value using Borland Turbo Pascal on an IBM AT compatible computer.
  160. The equation is due to Carl Friedrich Gauss (1777-1855).
  161.  
  162. In 1989, Pi was computed to 1,011,196,691 decimal digits.  This was done
  163. by the Chudnovsky brothers (David and Gregory), Columbia University,
  164. using super computers.
  165.  
  166. There has been several major brakethroughs that has allowed the
  167. calculation of Pi to get into the million digit range on a personal
  168. computer. A major limitation in Pascal is that array size is limited to
  169. 64K bytes and arrays use 16-bit integer as their index. In Borland C and
  170. C++, arrays are only limited by the amount of memory in your system and
  171. arrays can use long integer (0 - 2,147,483,647) as their index. This
  172. Pascal limitation forced me to convert my Multiple-precision packages
  173. from Pascal to C++.
  174.  
  175. The equations used to compute Pi were linearly convergent. This means if
  176. we want twice the accuracy, we need to compute twice as many terms at
  177. twice the precision. This was solved by using algorithm a or b.
  178.  
  179. The multiply, divide, and square root routines that would be used in
  180. algorithm a and b were running in time of the order of n squared O(n2).
  181. This meant that if it took 20 seconds to compute pi to 1000 digits it
  182. might take 230 days to compute it to 1,000,000 places. One way to change
  183. the order of the running time of the multiply is to use a fast Fourier
  184. transform (FFT) to perform the convolution of the digits between the two
  185. numbers and to perform the carry out of the digits after the
  186. convolution. This makes the running time of the multiply approximately
  187. O(n Log n). The divide routine can be changed to an iterative method
  188. that uses the multiply. The square root routine can be changed to an
  189. iterative method that uses the divide.
  190.  
  191. The last break through was getting an operating system that would allow
  192. you to easily use large arrays in extended memory. This was accomplished
  193. by using MS Windows and Borland C++ for Windows.
  194.  
  195. We now have a package on a PC that computes Pi in the 1,000,000 decimal
  196. digit range and runs in a time of approximately O(n Log3n).
  197.